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Q-t We present an efficient numerical method for evaluating Landau coefficients, which 

describe small-amplitude equilibrium states in the vicinity of the linear stability 
threshold. The method differs from the standard approach by the application of 
solvability condition to the discretized rather than continuous problem. Thus we 
avoid both the solution of the adjoint problem and the subsequent evaluation of the 
integrals defining the inner products in the standard approach. Instead of the ad- 
joint eigenfunction we use the left eigenvector of the discretized problem. The latter 
is supplied by the linear stability analysis together with the right eigenvector for the 
critical perturbation. Expanding equilibrium solution in small perturbation ampli- 
tude in the vicinity of the linear stability threshold, we obtain a matrix eigenvalue 
perturbation problem. Solvability of this problem requires its inhomogeneous term 
to be orthogonal to the left eigenvector. The method is demonstrated by using a 
Chebyshev collocation method to reproduce Landau coefficients for plane Poiseuille 
flow. 
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I. INTRODUCTION 



A number of fluid flows can become turbulent while being linearly stable. Some of these 
flows, such as, for example, plane Couette flow and circular pipe (Hagen-Poiseuille) flow, are 
linearly stable at all velocities, while some other may become linearly unstable at higher ve- 
locities. A typical example of the latter class of flows, which will be our main concern in this 
study, is plane Poiseuille flow. Theoretically, this flow is known to be linearly stable up to the 
critical Reynolds number Re c = 5722.22,- however, experimentally it has been observed to 
become turbulent at Reynolds numbers as low as 1000.-"- Moreover, turbulence is observed 
to develop in this flow on the time scales which are by several orders of magnitude shorter 
than those predicted by the linear stability analysis. 5 Such a subcritical instability can be ac- 
counted for by the positive feedback of perturbation on its growth rate, which is a non-linear 
effect beyond the scope of the linear stability theory. Thus, a perturbation of sufficiently 
large amplitude can acquire a large positive growth rate at subcritical Reynolds numbers, 
where all small-amplitude perturbations are linearly stable. For small-amplitude perturba- 
tions in the vicinity of linear stability threshold, this type of phenomenon is described, in 
general, by the so-called Landau equation. 6,7 Whether instability is sub- or supercritical is 
determined by the coefficients of this equation, which are referred to as Landau coefficients 
and have to be determined for each particular case. 

Basic formalism for deriving Landau coefficients for plane Poiseuille flow has been de- 
veloped by Stuart^ and Watson. 9 A more up-to-date account of this approach is given by 
Schmid and Henningson.— This method has been extended and modified by Reynolds and 
Potter— who used it to prove that plane Poiseuille flow is indeed subcritically unstable. Both 
aforementioned methods were compared by Sen and Venkateswarlu 1 ^ who applied them to 
calculate higher-order Landau coefficients for plane Poiseuille flow driven by a fixed pressure 
gradient. Not much difference between the methods was detected in the supercritical region, 
however the Watson method was found difficult to apply in the subcritical region, which 
represents the main interest for this flow. These type of asymptotic expansion methods 
have been reconsidered and surveyed by Herbert,— and substantially extended by Stew- 
artson and Stuart 14 who included a slow spatial variation, which resulted in the complex 
Ginzburg-Landau equation.— 

The evaluation of the Landau coefficients required in weakly non-linear stability analy- 
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sis is technically rather complicated. This may explain why most hydrodynamic stability 
stability problems are restricted to the linear analysis, which is of a limited practical signif- 
icance when the instability happens to be subcritical. A significant technical hindrance to 
the implementation of weakly non-linear stability analysis is the adjoint eigenfunction that 
needs to be found by solving the adjoint problem. Then several complex inner product in- 
tegrals containing the adjoint eigenfunction need to be evaluated in order to obtain Landau 
coefficients. 

In this paper, we develop a simpler but numerically more accurate method for evaluating 
Landau coefficients. The method is based on the application of solvability condition to the 
discretized rather than continuous problem. This allows us to evaluate Landau coefficients 
without using the adjoint eigenfunction, which in our approach is replaced by the left eigen- 
vector. Such a possibility has been briefly discussed by Crouch and Herbert.— A similar ap- 
proach employing Gaussian elimination has been noted earlier by Sen and Venkateswarlu.— 

The rest of the paper is organized as follows. In the section below, we formulate the 
problem for plane Poiseuille flow and introduce 2D traveling-wave solution. The latter is 
then expanded in small perturbation amplitude in the vicinity of linear stability threshold 
to re-derive the classic expressions for Landau coefficients. Section II III presents a detailed 
development of our approach for a Chebyshev collocation method, which is used to evaluate 
Landau coefficients for plane Poiseuille flow at its linear stability threshold. The paper is 
concluded by a summary of results in Sec. IIV1 

II. FORMULATION OF PROBLEM 

Consider a flow of incompressible liquid with density p and kinematic viscosity v driven 
by a constant pressure gradient Vpo = — g x Pq in the gap between two parallel walls located 
z = ±/i in Cartesian system of coordinates with the x and z axes directed streamwise and 
transverse to the walls, respectively. The velocity distribution v(r, t) is governed by the 
Navier-Stokes equation 



and subject to the incompressiblity constraint V ■ v = 0. Subsequently, all variables are 
non-dimensionalized by using h and h 2 /v as the length and time scales, respectively. In 
order to simplify the expressions of Landau coefficients obtained in the following, we employ 




(1) 
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the viscous diffusion speed u/h as the characteristic velocity instead of the commonly used 
maximum flow velocity. Due to this non-standard scaling Reynolds number appears as a 
factor at the convective term rather reciprocal at the viscous term. The problem admits a 
rectilinear base flow 



where u{z) = 1 — z 2 is the parabolic velocity profile and Re = XJ^hjv is the Reynolds 
number based on the centerline velocity Uq = 2Poh 2 /pu. At sufficiently high Re, the base 
flow can become unstable with respect to infinitesimal perturbations Vi(r, t), which due to 
the invariance of the base flow in both t and x = (x, y) can be sought as 



where A is the complex temporal growth rate for the Fourier mode with the wave vector 
k = (a, (3) and complex amplitude distribution vi(z). The linear stability analysis aims to 
find marginal values of Re depending on k for which neutrally stable perturbations defined 
by 9ft [A] = are possible. The lowest marginal value of Re, which is referred to as the 
critical Reynolds number, for this flow is Re c = 5772.22. It is attained at the critical wave 
vector k c = (1.02055,0), which corresponds to purely transverse perturbations defined by 
P = 0. For Re > Re c , the linear stability theory predicts exponentially growing perturba- 
tions. Evolution of these unstable perturbations depends on the nonlinear effects which may 
either inhibit or enhance the growth rate leading to what is known as super- and subcritical 
instabilities, respectively. The former is expected to set in only at supercritical Reynolds 
numbers, while the latter can be triggered by sufficiently large amplitude perturbations also 
in a certain range of supercritical Reynolds numbers. 

A. 2D equilibrium states 

In order to determine whether instability is super- or subcritical, we employ the approach 
of Reynolds and Potter, 11 known as the method of "false problems",—^ and search for an 
equilibrium solution in the vicinity of Re c as follows. The neutrally stable mode with 
a purely real frequency u> = — iA interacting with itself through the quadratically nonlinear 
term in Eq. ([T|) is expected to produce a steady streamwise-invariant perturbation of the 
mean flow as well as a second harmonic ~ e 2i(uit+ax) _ Subsequent nonlinear interactions 




(2) 




+ c.c, 



(3) 
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are expected to produce higher harmonics, which similarly to the fundamental and second 
harmonic travel with the same phase velocity c = —u/a. Thus, the solution can be sought 
in the form 

oo 

v(r,t) = E n v n (z), (4) 

n=— oo 

where E = e 1 ( ult + ax ) contains to, which, in general, needs to determined together v n by 
solving a non-linear eigenvalue problem. The reality of solution requires v_ n = v*, where 
the asterisk stands for the complex conjugate. The incompressiblity constraint applied to 
the nth velocity harmonic results in D ra • v ra = 0, where D n = e 2 J^ + ie x a n with a n = an 
stands for the spectral counterpart of the nabla operator. This constraint can be satisfied 
by expressing the streamwise velocity component 

u n = e x ■ w n = ia" 1 ?!/ (5) 

in terms of the transverse component w n = e z -v n , which we employ instead of the commonly 
used stream function. Henceforth the prime is used as a shorthand for d/dz. Note that Eq. 
(JSJ) is not applicable to the zeroth harmonic, for which it yields w = 0. Thus, u Q needs to 
be considered separately in our velocity-based approach. 

Taking the curl of Eq. (CQ) to eliminate the pressure gradient and then projecting it onto 



e y we obtain 



where 



and 



[D n - ium]Cn = K, (6) 



( n = e 9 -D n xv tt ={ (7) 

u' , n = 0. 



h n ^ ^ v?i— m ■ D m (^ m (8) 



are the y-components of the nth harmonic of the vorticity £ = V x v and that of the curl 
of the non-linear term h = V x (v ■ V)v. Henceforth, the omitted summation limits are 
assumed to be infinite. Separating the terms involving uq in sum (jHJ), it can be rewritten as 
h n = ia~ 1 (h'^ + h£), where 

K = n Yl ^{Wn-mVlw^ - w' m B 2 n _ m W n - m ) , (9) 

hi = ia n [u - UQB 2 n ]w n = M n {u )w n . (10) 



Eventually, using the expressions above, Eq. can be written as 

C n (iu, u )w n = h™, (11) 

with the operator 

C n (iu, uq) = [D 2 n - ium]D 2 n - Af n {u ). (12) 

The equation above governs all harmonics except the zeroth one, for which it implies w = 
in accordance with the incompressiblity constraint (151) . The zeroth velocity harmonic, which 
has only the streamwise component uq, is governed directly the x-component of Navier- 
Stokes equation (JT]): 

vZ = -Po + g Q , (13) 
where Pq is the dimensionless mean pressure gradient and 

90 = i a m™m™m ( 14 ) 

is the x-component of the zeroth harmonic of the nonlinear term g = (v • V)v. Velocity 
harmonics obey the usual no-slip and impermeability boundary conditions 

w n = w' n = Uq = at z = ±1. (15) 

B. Amplitude expansion 

The equations obtained above govern equilibrium states of 2D traveling waves of arbitrary 
amplitude. In the vicinity of the linear stability threshold, which represents the main interest 
here, the solution can be simplified by expanding in small perturbation amplitude. As 
discussed above, the basic harmonic ()3]) with a small amplitude 0(e) interacting with itself 
through the quadratically nonlinear term in Eq. (CQ) produces the zeroth harmonic, which 
modifies the mean flow, and the second harmonic, both with amplitude 0(e 2 ). These two 
harmonics further interacting with the basic one produce an 0(e 3 ) correction to the latter. 
The second harmonic interacting with the basic one also gives rise to the third harmonic 
with amplitude 0(e 3 ). This perturbation series is represented by the following expansion: 

oo 

w n = Yl e H+2 "M H \A\ 2m w nM+2m , (16) 

m=0 
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where eA = A is an unknown equilibrium amplitude of the basic harmonic and A = 0(1) is 
its normalized counterpart. The mean flow, which as mentioned above needs to be considered 
separately, is expanded respectively as 

u = u 0fi + e 2 \A\ 2 u ,2 + .... (17) 

Similarly, we expand also the Reynolds number and frequency 

Re = Re + e 2 Re 2 + ..., (18) 
tu = tu + e 2 Q 2 + ■ ■ ■ , (19) 

where Re$ is a marginal Reynolds number satisfying 5f[Ao] = for the mode W\ t x with the 
frequency uq = ^s[Xq] and the wavenumber a; e 2 Re2 = Re-2 and e 2 u)2 = u>2 are the deviations 
of Reynolds number and frequency from their linear stability threshold values. Substituting 
these expansions into Eqs. (TTT]) and (|T3|) . and collecting the terms at equal powers of e, we 
obtain the following equations. At 0(e°), we have the base flow equation 

<o = "^0,0, (20) 

where Po,o — 2i?eo and ^0,0 = i?eo(l— z 2 ) = Reou(z). At 0(e), we recover the Orr-Sommerfeld 
equation 

£i(iw ,M , )wi,i = 0, (21) 

which defines the linear stability threshold. Solution of this eigenvalue problem for a given 
wavenumber a yields Reo, Uq, W\^(z). The latter is defined up to an arbitrary factor, which 
for W\^(z) being an even function is determined by the standard normalization condition 

»i,i(0) = 1. (22) 

At 0(e 2 ), two equations are obtained 

< 2 = ~Po,2 - 2cr 1 SK* 1 U!y , (23) 
£ 2 (iw , u 0fl )w 2 ,2 = 2[(w 1 , 1 w' ljl )' - 2i&JJ', (24) 

which define the mean-flow perturbation 1*0,2 an d the second harmonic w 2 ,2 i n terms of 
Wi t i(z). The mean-flow perturbation depends also on the perturbation on the mean pressure 
gradient Po,2, which is zero when the flow is driven by a fixed pressure difference. Alterna- 
tively, if it is the flow rate rather than the pressure difference which is fixed, then P ,2 is an 



additional unknown which has be determined by using the flow-rate conservation condition 
/^i ^0.2 (z) dz = 0. We start with fixed mean pressure gradient Po,2 = 0. The the fixed flow 
rate case can readily be reduced to the former by including P ,2 into Re 2 as shown later on. 
To complete the solution, we need to proceed to the order 0(e 3 ), which yields 

A(iwo,Wo,o)rii,3 = + \A\- 2 [M l {Re 2 u + \A\ 2 u^ 2 ) + iu 2 D 2 ]w x , 1 , (25) 

where 

^i,3 = \ Ki D > 2 ,2 - w'^Blw^) - (w^Blw'*, - w'l^lw 2)2 ) . (26) 

Equation ( 125]) defines the correction of the basic harmonic in terms of the lower harmonic 
perturbations described above. It is important to notice that the l.h.s. operator of Eq. f l25|) 
is the same as that of the homogenous Eq. (12T|) . which is satisfied by W\^\. Thus, for Eq. 
( 125]) to be solvable, its r.h.s. cannot contain term proportional to W\^. Namely, the r.h.s has 
to be orthogonal to the adjoint eigenfunction : 

(u>t,i, Ks + \A\- 2 [Mx(Re 2 u + \A\ 2 u 0)2 ) + ic^D^^ = 0, (27) 

where the angle brackets denote the inner product. 10 This solvability can be rewritten as 

icu 2 = fnRe 2 + fx 2 \A\ 2 , (28) 

where 

fi 1 = -(wt )1 ,K(u)w 1 , 1 ), (29) 
^2 = -(^i,M(u ,2)^i,i + ^ 3 ) (30) 

for the adjoint eigenfunction normalized as (w^, D^u?]^ = 1. Equation (128]) represents a 
reduced Landau equation for the case of equilibrium solution, which requires u 2 to be real. 
This reality condition yields the sought equilibrium amplitude 

\A\ 2 = -Pe^H/^H, (31) 

which is the same as that resulting from the full Landau equation with the first Landau 
coefficient // 2 and the linear growth rate correction HiRe 2 .— Note that our non-standard 
choice of the characteristic velocity results in expressions (129]) and ( 130]) sharing operator J\f%, 
which simplifies their numerical evaluation in the next section. 
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The type of instability is determined by the sign of 9?[//2]- For an instability to be su- 
percritical, which supposes an equilibrium solution with |v4| 2 > to exist at positive linear 
growth rates Re2^[^i\ > 0, 9? [7/2] < is required. Otherwise, instability is subcritical. In or- 
der to calculate Landau coefficients and ( I30p following the standard approach outlined 
above one needs to solve not only the Orr-Sommerfeld equation (l2ip but also its adjoint 
problem for w^. Both the direct and adjoint problems, as well as those posed by Eqs. (1231) 
and (l24"j) . need to be solved numerically. Then the integrals in the inner products defining 
fii and fi2 need to be evaluated numerically. This standard approach can be significantly 
simplified by avoiding both the solution of the adjoint problem and the evolution of the inner 
product integrals. This is achieved by applying the solvability condition to the discretized 
problem as demonstrated in the following. 

III. NUMERICAL SOLUTION 

In this section, numerical evaluation of Landau coefficients will be demonstrated using a 
Chebyshev collocation method with Chebyshev-Lobatto nodes 

Zi = cos (m/N) , i = 0,---,N, (32) 

at which the discretized solution (w n ,Uo)(zi) = (w n ,Uo)i and its derivatives are sought. The 
latter are expressed in terms of the former by using the so-called differentiation matrices 
which for the first and second derivatives are denoted by Dfj and D^ 2 j. Explicit expressions 
of these matrices, which are too long to presented here, are given by Peyret.— Equations 
(I2ip . (123 p and <KZ%§ are approximated at the internal collocation points < i < N and 
the boundary conditions (I15jl are imposed at the boundary points i — 0, N. The operator 
£ n (iuj , m ,o) defined by Eq. (TT2]) . which appears in Eqs. (T2TT) and (1241) is represented by the 
matrix 

L n (iw , u ,o) = M n (u , ) - iwoA n , 

which contains 

M n (u 0l0 ) = F n [A 2 + Re N n (u)], (33) 
(A n )ij = (D^)j,j, 0<(i,j)<N, (34) 
where the latter represents the part of the collocation approximation of the operator 

(Pl)ij = Df] - al5 h] (35) 
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related with the internal nodes. The other matrix in Eq. (|33|) . 



(N n (u))jj = ia n [ui5ij - w"(A n )jj], (36) 

represents the collocation of approximation the operator defined by Eq. (1101) . Finally, the 
factor matrix— 

F n = I - B n (CA; 1 B n )- 1 CA- 1 (37) 

in Eq. (1331) is due to the no-slip boundary condition w'(±l) = 0, which is represented by 
Cw = with 

C l3 =Df], i = 0,N;0<j <N. (38) 
It also involves the part of the operator (135]) related with the boundary nodes: 

(B„) M - = (P 2 n )ij, 0<i<N,j = 0,N. (39) 

We start with the Orr-Sommerfeld equation, whose collocation approximation 

l_i(A, i2eu)wi,i = [Mi(i2eu) - AAi]wi,i = (40) 

after multiplication by A^ 1 reduces to the standard complex matrix eigenvalue problem 

[A^Mi^eu) - AIJwi,! = 0, (41) 

which can be solved using, for example, LAPACK's ZGEEV routine.— The marginal 
Reynolds number Re^ for a given wavenumber a is determined by the condition 9ft[Ao] = 
for the eigenvalue Ao with the largest real part. Simultaneously with the right eigenvector 
Wix, we find also the associated left eigenvector w^.— The right eigenvector is normal- 
ized using condition ( 122|) . and the left one is normalized against the former using the dot 
product of complex vectors wj 1 • w^i = 1. This normalization simplifies the expressions 
of Landau coefficients obtained below. Having found w^i, we can straightforwardly solve 
discretized counterparts of Eqs. (1231) and (l24p . which yield the mean-flow perturbation Uo,2 
and the complex amplitude distribution of the second harmonic w 2j 2- For the fixed flow rate 
considered later on, we need also the stream function of the mean-flow perturbation "00,2? 
which is obtained by solving the collocation approximation of ip' 2 = Mo,2 with the symmetry 
condition ^0,2(0) = 0. 
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Figure 1. Real and imaginary parts of the critical perturbation w\ t i given by the right eigenvector 
Wi l (a) and those of the respective left eigenvector wj x (b). 

Now we can proceed to solving the final equation (125]) . whose collocation approximation 
can be written similarly to Eq. (|40p as 

Lx(io;o, J Re u)w li3 = Fih^ 3 + \A\- 2 [F 1 N 1 (Re 2 U + |A| 2 u , 2 ) + ^Aijw^i, (42) 

which represents a matrix eigenvalue perturbation problem. For this system of linear equa- 
tion to be solvable, its r.h.s multiplied by A^ 1 as in Eq. (HTj) has to be orthogonal to w\ v — 
This discrete solvability condition leads to the same reduced Landau equation f l28|) . whose 
coefficients now are 

//! = -w t ljl -A- 1 F 1 N 1 (u)w ljl , (43) 
^ = -w^ • Ar 1 F 1 (N 1 (u , 2 )w lil + h? >3 ). (44) 

Owing to the symmetry of the problem, both w± t i and tto,2 are even, whereas 11)2,2 is an 
odd function of z. This allows us to search the solution in one half of the layer so reducing 
the number of required collocation points by half. M = N/2 = 32 collocation points in 
half layer is sufficient to obtain the the critical Reynolds number Re c = 5772.22, frequency 
bj c = —1555.18 and wavenumber a c = 1.02055 with six significant figures. 

The real and imaginary parts of the critical perturbation #1,1, which is given by the right 
eigenvector are shown in Fig. [T] together with the respective left eigenvector w{ l . Note 
that the latter is orthogonal to all other right eigenvectors but w^i and has only a numerical 
but no physical meaning. Because of different inner product definitions for the continuous 
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Figure 2. Velocity «o,2 an d the associated stream function ^0,2 of the mean-flow perturbation (a); 
the real and imaginary parts of the second harmonic amplitude W22 (b). 

and discrete problems, wj 1 is also distinct from the adjoint eigenfunction wf ± . Distributions 
of the mean-flow perturbation and that of the complex amplitude of the second harmonic 
are plotted for the top half of the layer in Fig. [2j Note that due to the non-standard scaling, 
our dimensionless frequency and velocity differ by a factor of Re c from the values obtained 
with the conventional scaling based on the centerline velocity. 

Substituting the above results into Eqs. (129]) and (130|) we obtain 

/ii = 0.0097118 - i0.222596, 
/i 2 = 0.0049382 -i0.0239131. 

As seen from Fig. El M > 32 collocation points produce Landau coefficients with about six 
significant figures. The first and most important result is 3ft [^2] > 0, which, as discussed 
above, confirms the subcritical nature of this instability in agreement with the previous 
studies. The coefficient /ii has been computed explicitly by Stewartson and Stuart 14 who 
found d\ = (0.17 + i0.8) x 10~ 5 for the standard normalization. Rescaling our result with the 
centerline velocity, we obtain fx\ = fii/Re c = (0.168251 — i3. 85633) x 10~ 5 , whose real part 
is close to that of di, while the imaginary is quite different. The reason for this difference is 
unclear. In addition, fix can be verified against the numerical results of the linear stability 
analysis for the complex growth rate in the vicinity of the linear stability threshold, where 
5X = A — A c ~ fii(Re — Re c ). As seen in Fig. HI the complex phase speed c = —iX/Rea, 
which is commonly used instead of A, is accurately reproduced by fii in the vicinity of Re c . 
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Figure 3. Relative variation of Landau coefficients with the number of collocation points M. 

In order to compare our Landau coefficient // 2 with the values found in the previous 
studies we have to take into account not only the non-standard scaling but also the fact that 
in our case A stands for the amplitude of the transverse velocity component w, whereas in 
the previous studies it denotes the amplitude of the stream function ip, which are related as 
w = icv0. Thus, our /i 2 rescales as 

jx 2 = ii 2 a 2 c Re c = 29.659 - H43.622. 

This result is close to /t 2 = ia c Ki = 29.46 — H43.41 found by Sen and Venkateswarlu^ 
using the method of Reynolds and Potter— for Re c = 5774, a c = 1.02 and c r = 0.2639. 
Note that K\ is mistaken for /i 2 by Schmid and Henningson^ who denote it by A 2 . Using 
a Chebyshev collocation method Fujimura^ found ^t[ji2\ = 30.962411, which appears to 
be less accurate than the result of Sen and Venkateswarlu^ obtained by a conventional 
finite-difference method. 

Reynolds and Potter— used their original method of "false solution" to obtain the first 
relatively accurate values of Landau coefficients for the case of fixed flow rate. Our solution 
for fixed pressure gradient can be converted to the fixed flow rate by using the non-zero 
pressure gradient correction Po,2 in Eq- (123]) . As seen from Eq. (1201) . this correction, which 
affects only the magnitude of the base flow, is equivalent to the substitution of i?e 2 by 

Re\ = Re 2 + \A\ 2 P 0t2 /2. 

The pressure correction Po,2, which according to the expression above produces a flow rate 
perturbation |A| 2 Po,2V>(l), has to compensate flow rate perturbation 2 1 ^4 1 2 -^o,2 (1) , which 
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Figure 4. Imaginary (a) and real (b) parts of the complex phase velocity c = — iX/Rea of the most 
unstable mode in the vicinity of the critical Reynolds number Re c calculated using m and supplied 
by the linear stability analysis (triangles) and taken from Sen and Venkateswarlu^ (dots). 

occurs at fixed pressure gradient. This results in 

Po, 2 /2 = -^ , 2 (1)M1) = 0.00217238, 

where t/>(1) = §q u(z) dz = |. Substitution of Re2 by Re\ in Eq. (128]) results in the replace- 
ment of fj,2 by 

= n 2 + = 0.0051492 - i0.0287487. 

Rescaling n\ with the critical Reynolds number based on the mean velocity Re c = \Re c = 
3848.08 and the critical wavenumber a c = 1.02071 used by Reynolds and Potter,— we obtain 

P% = v\ot\Re c = 20.64 - ill5.26, 

which is not far from fi\ = + ib^ = 19.7 — ill 1 obtained by Reynolds and Potter.— 

IV. SUMMARY AND CONCLUSION 

We have developed an efficient numerical method for evaluating Landau coefficients, 
which determine weakly non-linear evolution of small but finite amplitude perturbations in 
the vicinity of linear stability threshold. Our main interest was in the equilibrium states, 
particularly whether they exist above or below the threshold. This depends on the real part 
of the first Landau coefficient. If it is positive, the instability is subcritical. Otherwise, the 
instability is supercritical. 
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To identify the type of instability we employed the so-called method of "false solution" 
to search for equilibrium solution using expansion in amplitude, which was supposed to 
be small in the vicinity of linear stability threshold. This resulted in a matrix eigenvalue 
perturbation problem, which was solvable only when its inhomogeneous term was orthogonal 
to the left eigenvector. Application of the solvability condition to the discretized rather than 
the continuous problem was the main novelty of our approach. It allowed us to avoid both the 
solution of the adjoint problem and the subsequent evaluation of the integrals defining the 
inner products in the conventional approach. The adjoint eigenfunction in our approach was 
substituted by the left eigenvector, whereas the integrals were replaced by the dot products 
of complex vectors. The left eigenvector was supplied by the solution of the discretized linear 
stability problem together with the right eigenvector. 

The method was demonstrated by using a Chebyshev collocation method to reproduce 
Landau coefficients for plane Poiseuille flow at its linear stability threshold. Our results 
agreed well with the published values for both fixed-pressure-gradient and fixed-flow-rate 
driven flows. 

Circumventing solution of the adjoint problem not only considerably simplifies weakly 
non-linear stability analysis but also increases its accuracy by avoiding approximate evalua- 
tion of the inner product integrals. This makes our method more practical than the standard 
approach. 
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